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Abstract 



Simple application of the Einstein model combined with the elastic description of solid state is developed. 

The frequency of quantum oscillators has been assumed as volume dependent and, furthermore, elastic 

'■^ I energy terms of static character have been included to complete the description. Such an extension enables 

P3 ■ to construct the complete thermodynamics. In particular, the model yields practical equation of state 

O ' and describes the thermal expansion coefficient as well as the isothermal compressibility of solids. The 

^ v^^ , thermodynamic properties resulting from the Gibbs free-energy analysis have been calculated and illustrated 

in figures. Some comparison of the theoretical results with experimental data for solid argon has been made. 

►^ Keywords: Einstein model, equation of state, thermodynamic properties 

o\ _ 

^ ■ 1 Introduction 

(N 

.• ■ The Einstein model of the solid state has been known for many years as the first model able to describe qualita- 
r--s I tively the low-temperature behaviour of the specific heat [B [21 E] . This model often serves as an approximation 
fT^ ■ for studies of optical phonons, or soft modes in some intermetallic compounds H] and, as well, is found suit- 
able to describe thermal properties of some modern low-dimensional structures 5 . It is also referred to as 
a prototype for which the more sophisticated Debye model has been developed [5]. Nonetheless, it has been 
, ^ shown that in many cases the Einstein model provides better results than the Debye model [J. It is known 
rN ■ that in spite of its usefulness and simplicity, the pure Einstein model lacks the ability to describe fully the ther- 
J3 I modynamic properties of the solid state. For instance, neither the proper equation of state, nor the thermal 
" ■ ■ expansion or compressibility can be obtained. Thus, apart from the celebrated specific heat behaviour, which 
is an important consequence of the quantum nature of harmonic oscillators, other thermodynamic properties 
have not successfully been described within this model. 



There are numerous attempts in the literature aimed at improving the pure Einstein model [51 [91 [TOl [11] [12] . 
For instance, a "variational Einstein model" has been developed in Ref. [8] for describing low temperature 
solids from the Feynman path integral perspective. The formalism has been applied to a specific system, con- 
sisting of solid hydrogen with lithium impurities. However, one of the theoretical results of Ref. [8] seems to 
be controversial, namely the dependence of the free energy, which is presented there as an increasing func- 
tion of temperature (see Fig. 16 of Ref. [8 ), giving the indication that the corresponding entropy is negative. 
Cankurtaran and Askerov [9] have introduced the so called Einstein-Debye model for which the calculations 
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of the thermal expansion coefficient were possible by introducing the Griineisen parameter. However, neither 
the Gibbs free-energy, nor the compressibility could be calculated within that approach. Thus, the complete 
thermodynamic description has not been achieved there. In Ref. [TU] the "nonlocal Einstein crystal" has been 
considered, for which the thermodynamic functions have been constructed. The considerations presented in 
Ref. [TU] seem to be of pure theoretical character when the Einstein crystal is replaced by a harmonic sublattice 
with a size that is coupled to the total volume. Additionally, the Mie- Griineisen approach has been developed 
by Holzapfel et al. for the description of simple metals [TT as well as some simple solids, such as NaCl and 
MgO [12,. The optimized pseudo-Debye-Einstein model has been used there, and good agreement with the 
experimental data has been achieved. However, the method presented in Rcfs. [TTJ [T^] does not seem to be 
straightforward enough to use since it requires many fitting parameters to be introduced for calculations of the 
free-energy. 

On the other hand, there exists a rich literature on the description of solid state properties within the elastic 
models [131 [TU [TSl [HI [T71 [TS] [THl [201 [HI HI] • These studies are supported by the experimental measurements 
[531 [21] ■ For instance. Birch [T^] considers the single-crystal and polycrystalline NaCl at high pressure using 
the finite strain (Eulerian) theory. In this approach, the equation of state contains coefficients which have to be 
fitted in each particular temperature. Despite the usefulness of the model for studies of isothermal compression, 
it does not involve the thermal (vibrational) energy and therefore is not suitable for description of the specific 
heat. Similarly, in the paper of Vinet et al. [30] HT] the universal form of the equation of state has been found 
based on the scaling argumentation. The results obtained there have been discussed in the context of the 
Birch-Murnaghan equation of state [TB] and show very good agrement with the experimental data. However, 
the Gibbs free-energy has not been derived in that works (Ref. [20l \2T\ ). so that the full thermodynamic de- 
scription of the system has not been constructed. As a matter of fact, the same deficiency can be risen for some 
other papers devoted to the application of the equation of state based on the elastic models. A comparative 
review of the experimental equations of state existing in the literature is given in Refs. [^126) . 

Among other modern methods one should mention the first principles computations within the density- 
functional theory (DFT). This includes, for instance, the equation of state, elastic constants and phonon 
dispersion relation calculations [27]. On the other hand, by the first-principles molecular-dynamics simulation 
(FPMD) method the thermal expansivity and the specific heat have been calculated [5S]. A good agreement 
with the experimental data is often achieved. However, the analytic description remains always desirable from 
the point of view of understanding the physics behind the model. 

Such a situation motivated us to combine the Einstein and elastic models and to complete them in a way 
that enables to overcome the above-mentioned problems. Looking in the literature, such idea can already 
be found in some papers [301 130]; however, only particular thermodynamic characteristics have been studied 
there. Thus, our aim is to construct the self-consistent and full thermodynamic description. In order to make 
it possible, the following aspects should be taken into account: Firstly, the common frequency of quantum 
oscillators should not be treated as a constant value but should be volume (and, in consequence, temperature) 
dependent. Obviously, such idea is not new and for the first time emerged in the paper of Griineisen 'SP, where 
a specific frequency/volume dependence has been assumed as a hypothesis. This assumption is supported by 
the argument that the wavelength of collective excitations (phonons) depends on the crystal size. Secondly, 
the vibrational Einstein Hamiltonian necessarily needs to be completed by the static, elastic part [32]. The 
elastic energy accounts for the mutual interactions of the atoms even if they are not in a thermal vibrating 
state [IHl [22] . It is known that the static energy is responsible for crystal compressibility [T5] , and owing to 
its magnitude and nature, cannot be inferred from the model of independent oscillators. It turns out that for 
the purpose of the combined model the elastic free-energy should also be modified by taking into account the 
linear term (which is normally neglected in the elastic theory) as well as the static entropy. The details of this 
modification will be explained in the next section. Owing to analytical simplicity we will neglect the electronic 
excitation term. It is argued that this term is generally a small correction to the solid equation of state [30] . 

Taking this into account, we propose improving the model and constructing the Gibbs energy of the sys- 
tem from the beginning. We assume that the balance between the internal and external pressures keeps the 
system in mechanical equilibrium. According to our knowledge, the detailed balance between the expanding 
pressure of quantum oscillators and the compressive elastic pressure has not been discussed in the literature. 
Having obtained the expression for the Gibbs energy, a full thermodynamic description of the crystal can be 
achieved. The main parameters of the theory, which can be extracted from experiment, are the Einstein fre- 
quency in the ground state (or the Einstein temperature), the volume elastic modulus in the ground state - 
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supplemented by the structural space-filling coefficient - and the Griineisen parameter. We will show that these 
parameters, in principle, allow us to calculate other properties within a reasonable range of experimental values. 

The paper is organized as follows: In the following, theoretical section, the outline of the model is given 
and the method of derivation the Gibbs energy is presented. In particular, for the combined model the new 
equation of state is obtained. In the last section, some representative numerical results concerning the thermo- 
dynamic properties are illustrated in figures. These calculations are compared with the experimental data for 
solid argon. A critical discussion of the presented approach is also included. 



2 Theoretical model 

2.1 The model Hamiltonian 

The Hamiltonian is assumed in the form of: 

T~L = He + ^wj (1) 

where the elastic, volume dependent part, H^, can be written as: 

He^NAe+^NBe^-^^NCe^ + ^NDe^, (2) 



whereas the oscillatory part, T-Lun is given in the form: 



3Af 



i=l 



1 



n^^}_^hLulh, + -]. (3) 



In ^ and ([2]) iV is the number of atoms which are treated as the three-dimensional oscillators, and fii is the 
excitation number operator connected with the i-th oscillator. The form of the elastic Hamiltonian ^ assumed 
in this paper is different (much simpler) from that presented in Refs. [29] or [30]. The elastic Hamiltonian, He, 
consists of several terms. The most important is the harmonic term (ex e^j, where _B-constant is the volume 
elastic modulus in the ground state. We found that the linear term (oc e) is also necessary, where A-constant 
is responsible for the internal (compressive) pressure. Although the linear, anharmonic part is assumed to be 
very small in comparison to the harmonic one {A ^ B), its role is very important in balancing the internal 
(expanding) pressure produced by the anharmonic oscillators, i.e., by the 'Hcj part. A certain requirement for 
^-values is imposed as the equilibrium condition for the whole energy, which is discussed later. The relative 
elastic deformation s is defined by the relation: 

y = Fo (1 + e) , (4) 

where V is the crystal volume in a current thermodynamic state, and Vq is the volume at external pressure p = 
(vacuum) and temperature T — 0. Equation ^ indicates that a non-zero value of volume deformation e occurs 
when the external pressure p ^ and/or the temperature T > is applied. Thus, a possibility of introducing 
the isothermal compressibility and thermal expansion coefficient has been opened within the presented approach. 

A composition of the Hamiltonian ([T]) consisting of the classical {He) and quantum (T-Luj) parts is justified 
for the purpose of constructing the thermodynamics [22j . Obtaining the free energy of the full system is one of 
the main goals and that can be conveniently done by summing up the free energies corresponding to those two 
parts. A similar modus operandi has also been presented in the paper Ref. 32: albeit for a different approach. 
It should be stressed that by introducing the elastic Hamiltonian, He, the static no n- vibrational energy has 
been taken into account, which is crucial for the equation of state. It is worth mentioning that the method can 
be further generalized by introducing also the electronic part of the free energy i22j . 

2.2 The elastic free-energy 

In order to calculate the elastic (Helmholtz) free energy, one has to evaluate the static internal energy and 
the static entropy. The static internal energy Ue can be immediately found from the classical part of the 
Hamiltonian: 

Ue = He (5) 
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From Eq.(l2|) it is seen that this energy is explicitly volume-dependent. On the other hand, the static entropy 
Se is connected with filling up the volume of the system in the absence of thermal movement. Let us denote 
by Vex the own volume occupied by N atoms in the crystal state. For a given crystallographic lattice we can 
introduce the space filling coefficient q = T4a:/Vb. For instance, for FCC and HCP structures and the model of 
hard spheres q w 0.74. Thus, inside each primitive cell g = pi is the probability of finding an atom, whereas 
1 — q ^ P2 is the probability of opposite event. The static entropy connected with the bimodal distribution 
Pi (where i = 1,2) can be calculated from the mean value of Inpi, using the general formula for the exact 
differential of entropy: [55] 

dS,^-NkBd{liip,). (6) 



Hence, in our case: 



Se = -NkB [qlnq+ (1 - q) In (1 - q)] - NksC 



(7) 



where c is a constant of integration. Since c is arbitrary, as a matter of convenience we can choose c = 0. 
This arbitrariness is due to the fact that only changes in entropy are experimentally measurable and have 
thermodynamic significance [331 134| . 

The entropy is additive with respect to the number of primitive cells N. For the purpose of this paper it is 
assumed in the first approximation that the space-filling coefficient q is volume independent and, in consequence, 
the corresponding entropy S^ is constant. Thus, the static entropy (for c — 0) results in the residual entropy of 
the crystal for T — ^ 0, which is due to the fact that the whole volume is not perfectly filled. Moreover, the im- 
portant role of S'e (being dependent on q) may manifest itself in the structural (1st order) phase transitions [22j . 



The elastic, Helmholtz free-energy can now be found from the formula: 

F, = U,- TS,. (8) 

Hence, we finally obtain: 

Fe = 

+ NkBT[qlnq+{l-q)lnil-q)]. (9) 

For the full thermodynamic description this energy should be completed by the vibrational, Einstein part. 



2 3! 4! 



2.3 The vibrational free-energy 

For the Hamiltonian ([3]) it is convenient to employ the canonical ensemble. The Helmholtz free-energy can be 
found from the formula: 

F^^-kBTlnZ^ (10) 

where the statistical sum Z^ can be calculated exactly as: 



T?' e 



-pn^ 



2 sinh ( —(ilioj 



-3N 



where /3 = 1//cb2^- Thus we obtain: 



F^ ^iNkBTl-a 



2 sinh ( —fUhcj 



(II) 



(12) 



On the other hand, the vibrational internal energy Uuj is given by the statistical mean value of the corresponding 
Hamiltonian: 

u. = iu.) = mn. ((..) + 1) ^ \n^-J^^. (13) 

Hence, the vibrational entropy S^ can be found from the relationship: 

JJ ~ F 



s^ 



T 

^N Tiuj 

2T tanh {llShio) 



3NkB In 



2 sinh ( —(3huj 



(14) 
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The specific heat at constant V is given by the well-known formula: 

- 37VfcB(|.)'smh-(|,) (15) 

where Q = fiuj/k-Q. 

Further, in order to improve the Einstein model it is assumed that lj is not a constant but volume dependent 

according to the Griineisen assumption [31j : 

^ « -^ (16) 

where 7 is the Griineisen constant. Thus, we can write: 

Wo 



{l + ef 



(17) 



where we made use of the relation (|4]), and wq in (I17p is the Einstein frequency in the ground state (i.e., at 
p = 0, T = 0, and e = 0). Then, the variable 8 = huj/k-Q can be presented in the form of: 

where 

eo = ^. (19) 

00 is the so-called Einstein temperature. The correction of the Einstein model presented by Ea. p5|) implies 
that all the thermodynamic potentials will depend on the volume via elastic deformation e. 

2.4 Thermodynamic functions of the combined model 

The total Hclmholtz free-energy is given by the sum of expressions ([9]) and (fT2)) . with the help of (fT6|) - ([T9| : 

F = NAe+-NBe^ --NCe^ + —NDe'^ 
2 3! 4! 



A^fcBT[gln(7+(l-(7)ln(l-(7)] 
60 1 



+ 3A^fcBrin{2sinh 
The total Gibbs energy is then given by: 



2r(i + £)^ 



}• (20) 



Gip,T)=F+pV (21) 

The equation of state can be obtained from the variational principle: 

which, from (PTj) . is equivalent to 

nT/„ = - I 

de / J, 

where F is given by (OU)) . From (P^ . after differentiation of (pn|) . we obtain 



pVo = (^) (23) 



pVo = -NA-NBe+^NCe^ -^^NDe^ 

+ h — .J'^^ . , ■ (24) 



3 NkBQ 





2 ' (1 + ey^^ tanh 


[eo 1 1 

2T (l+e)T 



For the high-temperature limit equation (j24p takes a form: 



dF, 3NkBT 
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which is the classical lattice dynamics pressure [22]. On the other hand, in the low-temperature limit {T — > 0) 
and for small pressure the deformation e is small too, and Ea. (l24p can be linearized as: 



pVo 



-NA - NBe + --/NkBBo [1 - (7 + 1) e] . 



(26) 



Naturally, when p = then in the ground state is e = as well. Consequently, we have to demand that 
the internal pressure of oscillators, which then is equal to ^'^NkB&o, cancels out the pressure arising from 
the anharmonic, linear term, i.e., —NA. Then the system is in equilibrium without any external force. The 
condition for this demand leads to the specific choice for the anharmonic constant A, namely: 



where 



A = 37 Ao 



Aq = -hojo = -^fcBBo. 



(27) 



(28) 



The energy constant Aq presents a convenient unit for introduction of the dimensionless external pressure tt 
defined by: 

- ^° -P (29) 



and the dimensionless temperature t: 



r = 



Ao ■ 



(30) 



With the above quantities, and Eas. ipT)) - (j28p . the equation of state (j24p takes the simple form: 



TT -I- 37-1- -—e 
Ao 



1_C 
2Ao 

37 



1 D 

3! a;' 



(i+e) 



7+1 



tanh 



■(l + e) 



(31) 



Eg. (1311) . obtained within the approximations assumed in this paper, presents a dimensionless equation of state 
from which the elastic deformation e{7r,T) can be calculated for arbitrary temperature r and external pressure 

TT. 

The Gibbs free-energy (|2ip in the dimensionless form can be written as follows: 



G{p,T) 

NAn 



1 B o 1 C o ID, 



+ r[(7lng+(l-(?)ln(l-(?)] 
1 



3r ln{2sinh 



r(l + e) 



} + 7r(l + e), 



(32) 



where e = e (tt, t) is given by (j3ip . It is easy to show that from the formula p2p all thermodynamic properties 
can be derived self-consistently. For instance, the total entropy fulfills the relationship: 



dG\ 



fee fdG 



(33) 



and, with the use of pip , it can be obtained in the dimensionless form: 

S 



NkB 



-qlnq- (1 -g)ln(l - q) 



Til + ey 
31n{2sinh 



tanh ^ 

1 

[^(1 + 




1 




L + e)^J 
}■ 



On the other hand, the volume of the sample satisfies the equation 

V=[—] ^0 [dG 

\dpjj, NAo V dn 



(34) 



(35) 
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which again leads to the equation of state ([3T|). Other thermodynamic properties follow from the second 
derivatives of the Gibbs energy. For instance, the heat capacity at constant p is given by: 



Cn 



^-T(m .-- 



\9Ty, 



Ao 



9r2 



whereas the heat capacity at constant V can be calculated from the relationship: 



Cv 



^-Tf—] ^(—\ =^(^ 



\dTy, 



dT 



^0 V^^ 



(36) 



(37) 



On the other hand, the internal energy, [/, can be written in the dimensionless form: 



U 

'NAn 



376 - 



1 B 

2'M 



1 C 



1 D 



{l + ef 



tanh 



3! ^0 4! Ao 



1 



■{l + ef 



Application of (p8)) in (137)) leads to the dimensionless expression for Cy: 

2 



Cy 

NkT, 



1 



■(! + £)' 



sinh 



1 



■{l + ef 



(38) 



(39) 



It should be noted that for the particular case, when e = 0, Ea. ((M)) is equivalent to formula ([T5|) . However, in 
general, a non-zero e should be determined from the equation of state ([31]) . 

Among other response functions which can be calculated either from the Gibbs free energy, or from the equation 
of state is, for example, the thermal expansion coefficient: 



1 

V 



and the isothermal compressibility: 



kt = 

V 



dV_ 
dT 



dp / T 



fcs _ 

Ao 1 



Vn 



1 



NAo 1 






(40) 



(41) 



In particular, for p — and T ~ 0, i.e., in the ground state, these response functions can be obtained from the 
approximate equation of state (j26p . For this purpose. Eg. (1261) can be re- written in the form: 



Vn 



1 



■P- 



NAo3^{j + l) + B/Ao'^- ^^^^ 

It can be observed that Eq. (|l^ for p < presents Hookc's law at T = 0. Further, we can formally make use 
of the linear expansion of volume, namely: 



V{p,T)^Vo+(%) T + 

It can be noticed that Ea. (H5)) can equivalently be written as: 



dV\ 



P- 



(43) 



(44) 

where ao and kq are the response functions (I40p and (|4ip . respectively, and these functions are taken in the 
ground state. Now, by comparison of Eas. (|12|) and (jH]) we obtain the result: 



ao 



def J_ 

Vn 



and 



def 
Ko = 



1 

Vo 



dV\ 
Vo 







&v\ 

dp J (p=o) ~ NAo 37 (7 + 1) + B/Ao ■ 



(45) 



(46) 



It is apparent that the both elastic coefficients, A{= 37j4o) and B, are important for determination of isothermal 
compressibility at T = 0. Since B ^ Aq, where Aq is the zero-point energy of Einstein oscillators, it is obvious 
that for this compressibility an elastic energy plays a dominant role. Finally, let us note that the equation of 
state derived by Vinet et al. [H] for T = reduces for small p to the form: e « {l/Bo)p, which agrees with 
our eauation (H^ in the limit Ao -^ 0. 
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2.5 A note on the Griineisen equation 

The Griineisen parameter 7 is introduced by the definition [51 13] : 

which is equivalent to our relation (J16p . The typical experimental values of the Griineisen parameter belong 
to the range 1 < 7 < 3 [5]. It can be mentioned that some particular definitions of that parameter, depending 
on the model, have been discussed in Ref. 35 . The experimental determination of 7 is connected with the 
so-called Griineisen equation: 

-i^ = 7 -— (48) 

kt V 

which has been tested extensively, especially in the low temperature region. It is worth noticing that the 
derivation of (^5]) can be based on the exact thermodynamic identity: 

3fp (dp 



.r vary/ (^^) 

In the case in hand, in order to calculate {dp/dT)y , the equation of state p = p{T, V) should be taken in the 
form of ([M]) . It can easily be checked that after the calculation of the derivative in (P^)) . eq.(jl5]) is satisfied. 
It should also be noted that the calculations of all quantities in pS)) using the present method require prior 
solution of the equation of state (j3ip . In particular, for the low-temperature region, where Cy tends to zero 
according to the 3rd law of thermodynamics, we obtain from ea. (|15)) that ap -^ 0. This limit is in agreement 
with the previous result (I45[) . 



3 Numerical results and discussion 

3.1 Exemplary calculations for the model 

In order to perform the numerical calculations, based on the formulas from the preceding section, it is necessary 
to estimate the energy constants, from which ^0 and B are the most important. The volume elastic modulus 
in the ground state, B, can be found, for instance, from the sound velocity measurements and its experimental 
value is of the order (10^^^ ^ 10^^^) J. In turn, Ao-coefficient is given by Eg. (1551) . where Oq is the Einstein tem- 
perature. Assuming that 9o is typically of the order ~ (10^ -^ 10^) K we can estimate ^0 as ^ (10~^^ -;- 10"^") 
J. Thus, the A-constant {A = 37^0) is about 2 orders smaller than the elastic bulk modulus B, and the 
linear anharmonic term in the Hamiltonian ([2]) can be considered as a small correction to the harmonic one. 
This confirms the fact that such quantities as the thermal expansion coefficient ap (Eq. HO)) or the isothermal 
compressibility kt (Eq. I4ip are mainly determined by _B-constant, not by ^o- The above estimations of A 
and B constants, together with the assumed C = D = 0, yield from Eg. (HO)) ap ^ (10^^ -^ 10^"*) 1/K for the 
temperatures near the Einstein temperature, which is a physically reasonable order of value. In the same token 
assuming a realistic amount of volume per atom, i.e., Vq/N ~ (10~^^ -^ 10"^*) m'^, we obtain on the basis of 
Eg. ipi Kt ^ (10^^^ ^ 10^^") 1/Pa, which is also a correct order for the isothermal compressibility. 

Although the anharmonic elastic energy ex ^ is a small correction in comparison to the harmonic one, it can 
be of the same order of magnitude as the vibrational energy resulting from Einstein oscillators. Therefore, this 
anharmonic term is vital in balancing the expanding pressure resulting from oscillators, in the case when the 
frequency is volume dependent p7p . The expanding pressure arises owing to the decrease of the energy (fre- 
quency) of oscillators when they experience collective excitations in increasing volume. The total equilibrium of 
the system requires balancing of the internal (stretching) forces resulting from expanding quantum oscillators, 
internal compressive forces of the elastic medium and the external pressure p. The resulting deformation e is 
temperature and pressure dependent and can be calculated from the equation of state (1511) . The other energy 
constants, C and D, play merely a correcting role in modeling the static potential, and are important for high 
temperatures where the elastic deformation e is significant. The exemplary numerical results, presented in this 
section in Figs. 1-8, are obtained for the following set of parameters: B/Aq = 10^, D ~ and q = 0.74. 

In Fig.l we present the dependence of s upon tt for different reduced temperatures t — 0, 0.5 and 1, as 
well as for two selected Griineisen parameters: 7=1 and 7 = 2. The calculations are based on the equation 
of state (|3T|) for C/Aq — D/Aq = 0. It is seen that these dependencies are almost linear in character, and 
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Figure 1: Elastic deformation e vs. dimensionless pressure tt. The curves are selected for three dimensionless 
temperatures: r = 0, 0.5 and 1 and two Griineisen parameters 7 = 1 (dashed) and 7 = 2 (solid). The elastic 
energy parameters are : B/Aq = 100, C/Aq = D/Aq = 0. 
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Figure 2: Elastic deformation e vs. dimensionless temperature t. The curves are selected for three dimensionless 
pressures: n = 0, 0.5 and 1 and two Griineisen parameters 7=1 (dashed) and 7 = 2 (solid). The elastic energy 
parameters are : B/Aq = 100, C/Aq = D/Aq = 0. 

the slopes of the curves correspond to the isothermal compressibility. It can be deduced from Fig.l that the 
isothermal compressibility coefficient should be weakly dependent on pressure in a wide range of temperatures. 
This observation is in accordance with the analogous assumption in the Murnaghan theory of the equation 
of state p8]. At higher temperatures (for t — 1) the Griineisen parameter has remarkable influence on the 
magnitude of elastic deformation, leading to the increase of e when 7 increases. 



In Fig. 2 the dependence of e upon t for various external pressures (tt = 0, 0.5 and 1) is presented. The rest 
of parameters are the same as in Fig.l. In this case the dependencies are not linear and their local slopes are 
attributed to the thermal expansion coefficient. It is remarkable that by increasing the Griineisen parameter 7 
the relative deformation e increases, which is evidently pronounced at high temperatures. Such behaviour is in 
agreement with the previous figure (Fig.l). The dependence presented in this figure is qualitatively similar to 
the lattice constant dependence on temperature, as can be seen in the Ref. [36 for the case of AIN. The effect 
of external pressure on the temperature dependence of relative deformation in Fig. 2 resembles, for example, 
the results obtained in the Ref. [25] . 

In Fig. 3 the dimensionless thermal expansion coefficient (see Eal40p is plotted vs. temperature r for external 
pressure tt = and Griineisen parameter 7 = 2. Different curves correspond to various C/Aq parameters. It is 
seen that C/Aq has influence mainly at high temperatures, where it causes the increase of the thermal expansion 
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Figure 3: The dimensionless thermal expansion coefficient -r^ ap = j^ (^) vs. dimensionless temperature 
r. Different curves correspond to C/Aq = 0, 100, 200 and 300, whereas B/Aq = 100 and D/A^ = 0. The 
external pressure is tt = and Griineisen parameter is 7 = 2. 
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Figure 4: The dimensionless isothermal compressibility coefficient -^^^ kt = —jj- (f^)^- '^s. dimensionless 
temperature r. The parameters are the same as in Fig. 3. 



coefficient. Let us note, on the basis of Eq. ([25)1 and pop . that t — 2 corresponds to the Einstein temperature. 
In the low-temperature region the thermal expansion coefficient tends to zero independently on C/Aq, which is 
in agreement with the Griineisen equation (H51) . A qualitatively similar behaviour for the thermal expansivity 
has been found in Ref. fni for the case of MnO, in Ref. [TS] for Ti, Al, NaCl and Na, or in Ref. |3S] for AIN. 
The conclusion that the anharmonic term becomes important only in the high-temperature region is in agree- 
ment with the free energy calculations from first principles [25] . The similar conclusion can be drawn from the 
paper Ref. [3 7) where the influence of intrinsic anharmonicity on the thermodynamic functions has been studied. 



The dimensionless isothermal compressibility (see Eg HIT) is plotted in Fig. 4 vs. r for the external pressure 
TT — Q and Griineisen parameter 7 = 2. As in Fig. 3, different curves correspond to various C/Aq parameters. 
A remarkable influence of C/Ag-values on the compressibility is seen for the temperatures r > 1/2. On the 
other hand, for t — ^ the compressibility tends to some non-zero value, kq, being in agreement with Eg. ipS)) . 
A qualitatively similar experimental results have been obtained for Pb and NaCl. [18] The dependence of 
isothermal compressibility on temperature corresponds to the fact that the elastic constants are temperature 
dependent (the exemplary studies can be found in Refs. [5^IB5] V 



In Fig. 5 the Gibbs free-energy per 1 atom is presented in Ao-units vs. reduced temperature r. The three 
curves in Fig. 5 correspond to different external pressures: tt = 0, 0.5, and 1, and are plotted by the solid, dashed 
and dashed-dotted lines, respectively. In this case C /Aq = and 7 = 2. We see that application of the external 
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Figure 5: The Gibbs free-energy per 1 lattice site in ^o units upon dimensionless temperature r. Three values 
of the dimensionless pressure tt — 0, 0.5 and 1 are assumed. The Griineisen parameter is 7 = 2 and the elastic 
energy parameters are: B/Aq = 100, C/Aq = D/Aq = 0. 
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Figure 6: Entropy per 1 lattice site in Boltzmann constant units upon dimensionless temperature r. The dashed 
and solid lines correspond to 7 = 1 and 7 = 2, respectively, whereas BjA^ = 100, C j Aq = D j Aq = and 
TT = 0. The dashed-dotted line is plotted under constraint e = 0, i.e., for the ideal Einstein model. 



pressure causes the increase of the Gibbs energy. It is demonstrated in Fig. 5 that the Gibbs energy is a concave 
function of temperature with a monotonously decrease vs. r. Such behaviour indicates thermodynamically 
stable solution when the entropy (defined by Eg 155)) is positive. Let us note that the initial slope of the Gibbs 
energy, at r — > corresponds to the residual entropy. 

The entropy is illustrated in Figs. 6 and 7 as a function of temperature r. In both figures C jA^ — 0. In 
Fig. 6 we plotted entropy per 1 lattice site in the Boltzmann constant units, for external pressure tt = 0. The 
solid and dashed curves, which are for 7 = 1 and 7 = 2, respectively, correspond to our model, where £ is a 
function of temperature and is calculated from the equation of state ([5T|) . Due to our choice of integration 
constant (c = in Eq.Q) we see that the residual (configurational) entropy is present, which does not depend 
on the Griineisen parameter 7. Moreover, one can notice that for T ^ the entropy change vs. temperature 
tend to zero (AS" -^ 0), which is in accordance with the third law of thermodynamics. On the other hand, the 
increase of 7 causes some small increase of entropy at high temperatures. For comparison, the dashed-dotted 
curve presents the entropy for the pure Einstein model, where e is put equal to for all temperatures. As stated 
before, the condition e = is equivalent to the assumption that the frequency of oscillators is kept constant 
and does not depend on volume. For the pure Einstein model there is no residual entropy in the ground state, 
since the configurational ordering of atoms is not taken into account. 
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Figure 7: Entropy for the dimensionless pressure tt (tt = 0.5 and 1) normalized to the entropy for tt = 0, vs. 
dimensionless temperature r. The dashed and solid lines correspond to 7 = 1 and 7 = 2, respectively. The 
elastic energy parameters are the same as in Fig. 6. 
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Figure 8: Specific heat Cy at constant V per 1 lattice site in fee units vs. dimensionless temperature r. The 
elastic energy parameters are the same as in Figs. 6 and 7. In the main plot the parameters t: — Q and 7 = 2 
are assumed. In the inset the difference Cy (tt) — Cv (e = 0) is plotted for tt = and tt = 1 , as well as for 7=1 
(dashed line) and 7 = 2 (solid line). 
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It turned out that the entropy is very weakly sensitive to the external pressure, and in Fig. 7 we want to 
study this property in more detail. The relative changes of the entropy S{tt)/S{tt = 0) are plotted vs. r for 
TV = 0.5 and 1. The dashed and solid lines are for 7 = 1 and 7 = 2, respectively. We note that the sensitivity of 
S" on TT is greater for higher Griineisen parameter 7. It is also seen that the influence of high external pressure on 
the ratio S{n)/S{n = 0) is most remarkable for some intermediate temperatures (r w 0.6), whereas for r ^^ 
the influence of tt becomes negligible. Again, the behaviour of 5* near T = is in agreement with the third 
law of thermodynamics. The lowering of entropy for T > 0, when the external pressure is applied, corresponds 
to the positive thermal expansion coefficient, in accordance with the Maxwell relation: {dS/dp)j, — — [dV/dT) . 

The specific heat per 1 lattice site in the Boltzmann constant units and at constant volume is plotted vs. 
temperature t in Fig. 8. The anharmonic parameters are C/Aq — D/Aq = 0. The main plot has been obtained 
for TT = and 7 = 2. We note that the value of t = 2 corresponds to the Einstein temperature. The low- 
temperature behaviour of the specific heat is typical for the Einstein oscillators, being in agreement with the 
3rd law of thermodynamics. On the other hand, the high-temperature part of the curve is classical. In order 
to see the difference between the present model and the classical Einstein result in more detail, in the inset, 
we plot this difference for tt = and tt = 1, as well as for two Griineisen parameters: 7 = 1 (dashed line) and 
7 = 2 (solid). Cv'(7r) corresponds to the specific heat of the present model, whereas Cy(e = 0) is the Einstein 
result, when w does not depend on F. It is worth noticing that for tt = some small increase of the specific 
heat occurs for temperatures T > 0. Similarly to the entropy, the specific heat is only weakly sensitive to the 
external pressure tt. In particular, one can see that the pressure of tt = 1 causes a small decrease of Cy at some 
restricted range of low temperatures. When the Griineisen parameter increases all the above changes become 
enhanced. 

3.2 A comparison with experimental data for solid argon 

For solid argon, which forms FCC structure below the melting temperature of 84 K, the Debye tempera- 

1 /^ 

ture is 0d=85 K. From the approximate formula 0o/6d = (^r/B) ' , which can be derived on the basis of 
Ref. [5], the Einstein temperature, 9o=68.51 K, can be estimated. Hence, on the basis of eg. ipS)) we find that 
the Ap-constant amounts to Aq = 0.04729 x 10~^° J. On the other hand, the isothermal compressibility at 
zero temperature, kq, can be found from the experiment |39) . and for polycrystalline samples it amounts to 
kq = 4.1 X 10~^*'Pa~^. The volume per atom can be estimated as Vq/N = 36.383 x 10~^°m^. The Griineisen 
parameter can be taken as a mean value of the experimental data from various temperature ranges, |39j which 
yields 7 = 2.5. Having obtained the data given above, on the basis of eq. pS]) the ratio B/Aq can be estimated, 
yielding B/Aq — 160. The other elastic energy constants, which are of higher order, i.e., C and D, can be 
treated as the theoretical fitting parameters. We found that the optimal fit to the experimental data (taken 
from table HI, page 561 of Ref. [^) is obtained for C/Ao=1250 and D/Ao=5000, valid for e > 0. Taking into 
account the above set of parameters, we have calculated e, kt, cup and Cy for solid argon in the full temperature 
range < T < 84 K. To complete the comparison, we supplemented the results with two example isotherms. 
Since the isotherms cover the range of both positive and negative deformation e, we decided to allow some 
asymmetry in our model elastic potential. Thus, for e < we adopted the anharmonic parameter value of 
C/Ao=3000, i.e. the potential is more steep in this range. The results are presented in Figs. 9-13 by the solid 
lines, whereas the experimental data [39] are shown by the open symbols. In addition, the figures have been 
supplemented with the smoothed experimental data from the Refs. [lOlBTlH^ . 

In Fig. 9 the calculation of elastic deformation e vs. absolute temperature T is presented. The derivative 
of elastic deformation over the pressure, i.e., the isothermal compressibility coefficient, kt, is plotted vs. T 
in Fig. 10. It can be noted that in both figures (9 and 10) the experimental dependencies are nonlinear, in 
agreement with the present theory. In turn, the thermal expansion coefficient, a^, at constant pressure (p = 0) 
is plotted in Fig. 11 for the same temperature range as in Figs. 9 and 10. Finally, the molar specific heat Cy 
at constant V is presented in Fig. 12. For comparison, in Fig. 12 the dashed line is plotted for the pure Einstein 
model when we impose the constraint e = 0. 

Our calculated isothermal compressibility, kt, can be related to the adiabatic compressibility, ks-, via the 
formula: 

KT = Ks (1 + apjT) . (50) 

It can be checked that the above formula is equivalent to the Griineisen equation (PSJ) which is also satisfied 
in our case. As pointed out in Ref. [39], the values of ks and kt can be measured independently; ks by 
the ultrasonic method from isentropic sound velocity, and kt by the piston-displacement method of Bridgman. 
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Figure 9: Elastic deformation e vs. absolute temperature T for solid argon. Solid line - our calculation, dashed 
line - smoothed experimental data after Ref. [ID], the open symbols - experimental data after Ref. [55] . 




Figure 10: Isothermal compressibility coefhcient kt vs. absolute temperature T for solid argon. Solid line 
our calculation, the open symbols - experimental data after Ref. |39) . 



o 




Figure 11: Volume thermal expansion coefficient Up at constant pressure p = vs. absolute temperature T 
for solid argon. Solid line - our calculation, dashed line - smoothed experimental data after Ref. [10], the open 
symbols - experimental data after Ref. 39 , the filled symbols - experimental data after Ref. [^ . 
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Figure 12: Molar specific heat Cy at constant V vs. absolute temperature T for solid argon. Solid line - our 
calculation, dashed line - smoothed experimental data after Ref. [ID], the open symbols - experimental data 
after Ref. [W. The dashed-dotted line correspond to the calculation with constraint e = 0, i.e., for the pure 
Einstein model. 




Figure 13: Isotherms for solid argon (dependence of elastic deformation on external pressure) for two selected 
temperatures. Solid lines - our calculation, dashed lines - smoothed experimental data after Ref. [JT]. 
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Using the formula (150)) a satisfactory agreement between those two measurements has been obtained in Ref. [39] . 
Such consistency directly transfers to our calculations, therefore in Fig. 10 only kt curve has been presented. 

Regarding the isothermal compressibility kq (at p = 0, T = 0) we adopted here the value kq = 4.1 x 
-|^Q-io pg^-i fQj. polycrystalline samples after Dobbs and Jones [39 . On this basis the classical (atomic volume) 
bulk modulus, B' , can be estimated at T = K: B' = I/kq = 0.244 x 10^° N/m^. On the other hand, for 
single cubic crystals the bulk modulus can be obtained from the formula 2 B' = (Cn + 2C12) /3, if the elastic 
constants Cn and C12 are known. These constants for argon single crystals have been deduced from measured 
isentropic sound velocity [13]. The values extrapolated to T = K are: Cn — 4.39 and C12 — 1.83 in units 
of 10^ N/m^ [43] . Hence, the bulk modulus B' obtained on this basis is B' = 0.268 x 10^° N/m^. It can be 
concluded that this value of B' for single crystal is about 10% larger than for polycrystalline samples given in 

In the Fig. 13 two selected isotherms are shown (low- and high-temperature one), presenting elastic defor- 
mation £ as a function of external pressure p. It is visible that using a model asymmetric potential, a good fit 
to experimental data in whole pressure range studied is obtained, both for high- and low-temperature results. 
This emphasizes the pronounced importance of the exact form of the static lattice potential for pressure studies, 
since the deformations here exceed 10%. The selected potential provides reliable description for pressures up 
to 10 kbar. In light of the assumption that our static elastic potential originates from the expansion at the 
point T — Q and p = 0, further adjusting the potential would improve the consistency between calculations and 
experimental data for higher pressures. 

It can be concluded, on the basis of Figs. 9-13, that the agreement between calculated lines and the 
experimental data is quite satisfactory. One should take into account that all these curves have been calculated 
selfconsistently for the same set of parameters, as described above in this subsection. It can be supposed that the 
fitting could even be better if we would allow the Griineisen parameter to vary with temperature and pressure, 
as it has been suggested from the experimental measurements [321111]. As far as the specific heat is concerned, 
the new result (solid line) differs only insignificantly from the Einstein result (dashed), and the difference is 
mainly noteworthy at high temperatures. However, it should be noted that the dashed line in Fig. 12 is the only 
result of the pure Einstein model which can be compared with all these experimental data. The full description 
of temperature dependencies of remaining quantities, such as e, kt and ap has been possible in the improved 
approach to the Einstein model, when the equation of state pip is taken into account. 

3.3 Final remarks 

In the present paper a simple combination of the Einstein and elastic models of solid state is presented. First of 
all, the dependency of frequency of quantum oscillators on the volume has been introduced in a simplified way 
via Griineisen assumption. Simultaneously, the elastic properties of the crystal have been taken into account 
via classical Hamiltonian, containing anharmonic terms. 

An idea that the free-energy is a sum of several components has been presented, for instance, in the book of 
Wallace i_22 . However, we have shown that exploitation of this idea led us to the new form of the equation of 
state (|3T|) . Contrary to the Birch-Murnaghan equation of state, which presents only isothermal description [22j . 
in our equation the temperature T plays a role equivalent to the rest of variables, i.e., p and V . In particular, 
the thermal expansion coefficient, as well as the isothermal compressibility have simultaneously been derived 
from this new equation of state. It is well-known that these quantities could not be inferred from the sole 
Einstein model. On the other hand, our equation of state includes also the pressure resulting from expanding 
quantum oscillators. Moreover, comparing our method with the papers based on the Wallace approach, one 
should notice that one of our anharmonic parameters in the static potential, namely A, is not independent, 
but has been related to the Einstein temperature (via Eas. [29land l30t . As we pointed out, such relationship 
assures that the equilibrium condition for the total free energy at p = and T = is satisfied and the system 
of quantum oscillators becomes stable. When this point is not discussed, the functional form of the static 
lattice energy can also be adopted in the form given in Ref. [5D] (after Vinet et al.). That form is much more 
complicated than our polynomial approach, nevertheless it has successfully been applied for very high pressures. 

The present method requires several constant parameters, such as: the Einstein temperature, isothermal 
compressibility (or the volume elastic modulus) at the absolute zero temperature, as well as the Griineisen 
parameter, which should be taken from experiment. 

The region of applicability of the model (its temperature and pressure range of validity) is mainly connected 
with the number of terms which are taken into account in the polynomial form of elastic Hamiltonian ([2]) . This 
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form is based on the assumption that the equilibrium central point is e = for T = and p = 0. Thus, the 
theory is best applicable for the expansion around this equilibrium. However, as we have seen on the example 
of solid argon, merely first four terms were sufficient to describe deformation e in a wide range; starting from 
e pa 0.1 (near the melting temperature and p = 0, Fig. 9) down to e w —0.16 (for T = 4 K and p = 10 kbar. 
Fig. 13). Of course, validity of the model for lower values of e < (for very high pressure, where anharmonic- 
ity of static potential plays a role) would require higher order terms in ([2]). One should also remember that 
our assumption concerning the space filling coefficient (q = const.) would not be fulfilled at extremely high 
pressures. As far as the vibrational anharmonic effects are considered, they all are taken into account by the 
effective Griineisen parameter 7. We think that the approach is useful in the range of temperatures up to the 
melting point. In conclusion, the presented formulation allows for a complete thermodynamic description of 
the system in a relatively wide range of external pressure and temperature. Our considerations are related to 
the quasistatic processes. Therefore, the shock-wave experiments cannot be described within this model. 

It should be noted that the method is relatively simple, gives analytical form of the equation of state, and 
therefore it can serve as a first approximation for more advanced approaches. For instance, in the prototype 
calculations we have used only single variable e for description of the volume elastic deformation. However, 
it seems possible to generalize the approach for anisotropic deformations and anisotropic external pressures, 
involving also the Poisson coefficient. The presented method can be potentially extended basing on the Debye 
approximation, in which the linear dispersion relation for collective excitations together with the proper density 
of states are taken into consideration. For instance, in Ref. |30) the vibrational free energy has been taken into 
account in the high temperature Debye model. In that approach each moment of the density of states function 
requires a separate Griineisen parameter, which makes the method much more complicated. On the other 
hand, for metallic systems the electronic part of the free-energy should also be included [30]. However, such 
extensions of the method need further studies and should be a subject of separate assignment. 

The paper has been partly supported by the grant VEGA 1/0431/10. 
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